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We re-examine the evidence of hemispherical power asymmetry, detected in the cosmic microwave back- 
ground (CMB) WMAP (Wilkinson Microwave Anisotropy Probe) data using a new method. We use a data 
filtering, preprocessing , and a statistical approach different from those used previously, and pursue an indepen- 
dent method of parameter estimation. First, we analyze the hemispherical variance ratios and compare these 
with simulated distributions. Secondly, working within a previously proposed CMB bipolar modulation model, 
we constrain the model parameters: the amplitude and the orientation of the modulation field as a function of 
various multipole bins. Finally, we select three ranges of multipoles leading to the most anomalous signals, and 
we process corresponding 100 Gaussian, random field (GRF) simulations, treated as observational data, to fur- 
ther test the statistical significance and robustness of the hemispherical power asymmetry. For our analysis we 
use the Internally-Linearly-Coadded (ILC) full sky map, and the KQ75 cut sky V channel foregrounds reduced 
map of the WMAP five year data (V5). We constrain the modulation parameters using a generic maximum a 
posteriori method. 

In particular, we find differences in hemispherical power distribution, which when described in terms of 
a model with bipolar modulation field, exclude the field amplitude value of the isotropic model A — at 
confidence level of - 99.5% ( - 99.4%) in the multipole range £ e [7, 19] (I e [7, 79]) in the V5 data, and at 
the confidence level ~ 99.9% in the multipole range £ G [7, 39] in the ILC5 data, with the best fit (modal PDF) 
values in these particular multipole ranges of A = 0.21 (A = 0.21) and A = 0.15 respectively. 

However, we also point out that similar or larger significances (in terms of rejecting the isotropic model), 
and large best-fit modulation amplitudes are obtained in GRF simulations as well, which reduces the overall 
significance of the CMB power asymmetry down to only about 94% (95%) in the V5 data, in the range £ £ 
[7, 19] (£e [7,79]). 



I. INTRODUCTION 



The Gaussianity and the statistical isotropy of the fluctuations in the Cosmic Microwave Backgrounds Radiation 
(CMBR) are two generic features of current standard cosmological model and are compatible with the simplest infla- 
tionary scena r ios. T hese predictions have been extensivel y studied in number of works. An inc o mplet e list includes: 
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2001h : lYadav & Wandeltl d2008l) and references therein. Within the theory of inflation the primordial fluctuations are expected 
to form a Gaussian Random Field (GRF) at the leading order in perturbation theory. Their statistical properties are imprinted 
in the CMB fluctuations, providing an interesting window on the processes of the early Universe. Although the instrumental 
effects, like non-Gaussian and non-isotropic noise, or eccentric beams, and astrophysical foregrounds effects, like Galactic, and 
extra-Galactic point sources, and extended sources of emission, are either well controlled, or corrected for, or mask ed out, a set of 
an unexpected anomalies of various magnitudes and at various scales have been detected in the current CMB data ( Abramo et all 
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Land & Magueiid 12005 allblld l2007h (see also ICabella & Marinuccil d2008l) : iHutererl d2006h : iMartinez-GonzalezI d2008l) for re- 



cent reviews and references therein). These anomalies call for plausible theoretical explanations since, if robustly detected, 
these c an be used as valuable observables of the physics of the early Universe, or a new window on some of the late time 
effects dAkofor et al l|2007l: IBernui & Hip61ito-Ricaldill2008l: iBrown & Crittendedl2005l: IDestri et al.ll2008l: iDvorkin et all2008l: 
lErickcek et al.ll2008l: llnoue & Silkll2007l) . 

In this paper we re-investigate the well-known, hemispherical power asymmetry observed in the CMB maps. We revisit the 
properties of the asymmetry, constrain parameters of the previously proposed bipolar modulation field model dGordon et al.l 
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120051) responsible for generation of the asymmetry, and we estimate the statistical significance using a generic maximum like- 
lihood method, and realistic Monte-Carlo simulations. Given that we introduce and utilize a different method, from those 
previously used, and rely on different assumptions, while relax some other, our results can serve as a separate cross consistency 
and stability check. We provide a through tests of the method so as to validate the presented results. 

For the first time we estimate the parameters of the hemispherical modulation for selected ranges of multipoles. Finally, we 
assess the significance of the hemispherical power asymmetry via direct comparison to the GRF CMB simulations. 

The main differences from, and extensions to the previous analyses are: (i) we rely on the local real-space measurements of 
the variance as an estimator for the power asymmetry, (ii) we do not assume any priors on the probability distribution function 
(PDF) for any of the modulation parameters, and explore the likelihood function in the full (albeit sparse) parameter space and 
apply interpolation, (iii) we fully include the effects of the cut-sky, cross-multipoles power leakage and (iv) we analyze the 
power asymmetry in the various slices through the spherical harmonic space, pre-filtering the data prior to the analysis, rather 
than considering all scales scrambled together. 

Since the full exploration of the parameter space is CPU expensive, our analysis is based on a few assumptions that greatly 
simplify, and speed-up the parameter estimation process. Using this different approach, while providing tests and justification 
for the assumptions made, we give a new estimates on the significance of the hemispherical power asymmetry anomaly. 

We also discuss limitations of usage of the method with regard to the extent to which the assumptions of the method remain 
acceptably valid. 

The organization of the paper is as follows: In Section |IJwe describe our datasets and CMB simulations. In Section [TTT1 
we present results of a statistics that measure the hemispherical power ratios. In Section [iVlwe focus on the properties of the 
power modulation model, our assumptions and tests of the assumptions, and then detail on our method for modulation parameter 
estimation. In Section [VI we present the results of various tests of the method. Results of the application of the method to the 
real CMB data are presented in section [VII Discussion and conclusions are given in sections Ivnl and rVTlTl respectively. 

II. DATA AND SIMULATIONS 

For the main analysis in the paper we use the WMAP five-year foreground reduced CMB temperature maps dHinshaw et al.l 
120081) from differential assemblies (DA) VI and V2, because these spectral channels provide the best trade off between fore- 
grounds of different spectral properties (i.e. the blue tilted galactic dust emission, and red-tilted galactic synchrotron and free-free 
emissions). We co-add these observations (and corresponding simulations) using the inverse noise pixel weighting scheme. We 
will refer to the resu lting map as V5. We generate simulations using the fiducial best fit AC DM model power spectrum of 
iDunkley et alJ d2008l) (constructed using the mean likelihood parameters) which we call Cf d . 

Also, for comparison purposes, we will use the five year release of the I nternally-Linear ly-Coadded (ILC5) map and also for 
additional tests the Harmonic-Internal-Linear-Combination (HILC5) map dKim et 

In Fig. [T] we plot the power spectra of the data sets that will be used in the power modulation parameters estimation analysis. 
For the purpose of the analysis (to be explained latter) we create a fitted power spectrum to the ILC5 data by concatenating 
the Cf d power spectrum in the limit of low multipoles {I < 30), with the cubic spline fit to the piece-wise averaged, full-sky 
power spectrum, reconstructed from the ILC5 data for multipoles i > 30. We will call this fitted spectrum C] LCfit . We cut 
off the residual strong foregrounds in the ILC5 map at the ±350/iK threshold; a level estimated from GRF foregrounds-free 
simulations. 

Note that in the analysis of the modulation parameter estimation we will focus only on the large scale multipoles (£ < 80), 
where the differences in the power spectrum due to the 350/iK chop of the residual galactic contamination remaining the ILC5 
data are completely unimportant, and where the signal-to-noise (SNR) ratio yields SNR > 10 2 (Fig.QJ. 

Throughout the paper we will use the V5 data along with the KQ75 sky mask unless otherwise noted. For comparison 
purposes, the ILC5 and HILC5 datasets are used without any masks throughout of the paper 1 . 



Although we realize that the usage of the unmasked ILC maps can lead to somewhat overestimated constraints on any power asymmetry signals due to residual 
foregrounds contamination, we include the analysis of these maps in the full considered multipole range (£ < 80) mostly for comparison purposes. 
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Figure 1: Power spectra of the data sets used in the analysis. We generate the GRF simulations of the V5 data using the best fit minimal 
ACDM model (green line) to the reconstructed WMAP5 data power spectrum (red line). The corresponding cut- sky pseudo power spectrum 
from one of the V5 simulations and V5 noise power spectrum are plotted in black (the top-black curve and the bottom-black curve respectively). 
Simulations of the ILC5 data are created using the fitted power spectrum C* LCfit (light blue). The reconstructed from the ILC5 data power 
spectrum is plotted using the dark blue (top) line. The noise power of the noisiest channel of the WMAP data (W) smoothed with one-degree 
(FWHM) Gaussian beam is also plotted (dark-blue bottom) and used here to place an upper limit constraint on the amount of noise in the ILC 
data. The power spectrum of the ILC5 data chopped at temperature threshold of 350/xK is also shown (pink line). In the parameter estimation 
analysis we will use only the range of multipoles £ < 80. 



These are simply the maximal (^NS,^ max ) and the minimal (TNS/ max ) ratios of the hemispherical standard deviations, found in 
the all sky search over directions n s , that define the orientation of two hemispheres. The on and as values are the cut-sky (in 
case of V5 data) or full- sky (in case of ILC5/HILC5 data) hemispherical standard deviations of a map. 

To define the grid of direction s (n s ), we choose to use the first 96, ring-ordered directions, defined by the pixel centres of the 
HEALPIX pixelization scheme dGorski et al.ll2005l) of resolution parameter n s = 4. The measurements are performed using 
V5, ILC5 and HILC5 datasets within either, chosen ranges of multipoles, or as a function of a cumulative maximal multipole 
number £ max . 

The results of this survey are summarized in Fig. [2] 

The hemispherical asymmetry in the filtered bins of multipoles seems to be localized within the range of multipoles 
4ssym e [8, 15] in all of the dataset: V5, ILC5 and HILC5. 

Note that the first point in every plot: i.e. the quadrupole (I = 2), has the ratios always equal unity, because only the 
quadrupole map is used in variance measurements, and since the single-multipole maps have a point (anty) symmetry, due to the 
properties of the spherical harmonics, the variance is identical in the two hemispheres. In the range of multipoles i G (2, 30) we 
process the maps containing only two neighboring multipoles: eg. for the £ = 3 we use combined maps of multipoles £ = 2 and 
£ = 3; for £ = 4 we use combined maps of multipoles £ = 3 and £ = 4, and so on. For higher multipoles (£ > 30) the bins are 
larger, and are defined by the two neighboring plotted points in Fig. [2 In the case of the cumulative plots (second and fourth 
row in Fig. [2]) we use all multipoles from £ = 2 up to £ m &x, and due to the cumulative process, the curves on the right-hand 
side plots, (showing the statistics in the large-^ limit), do not exhibit much of details, since most of the map power (variance) 
comes from the low multipoles. In the limit of large multipoles, the multipole range- filtered maps show much more details as 
they are not overwhelmed by the power of the low multipoles. In particular, notice the strong, systematical deviation away from 



III. HEMISPHERICAL POWER RATIO 



We begin the analysis of the hemispherical power asymmetry by computing the following statistics: 
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Figure 2: Hemispherical power asymmetry in filtered, in the spherical harmonic space, V5 data (top panel) and in the ILC5 and HILC5 data 
(bottom panel). In the first row of each panel we plot the maximal (i?NS,^ max ) (blue)/minimal (rNS,£ max ) (red) ratios of the standard deviations 
as a function of the filtered range of multipoles defined by the two £ values corresponding to any of the two neighboring plotted points (see 
text for details). In the second row of each panel we plot the ratios as a function of the cumulative maximal multipole number ^ max defining 
the considered range of filtered multipoles: £ G [2, £ ma x]. The red (dashed) line is the inverted minimal (l/^NS,€ max — r SN,£ max ) rat i° °f me 
hemispherical standard deviations, plotted to facilitate appreciation of the power asymmetry by direct comparison with the blue lines. The gray 
bands represent the 68%, 95% and 99% confidence level contours. The simulations' mean is plotted with yellow solid line. The hemispherical 
asymmetry mostly appears to be confined to the multipoles range 4ss ym G [8, 15] in V5, ILC5 and HILC5 dataset. 
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the simulation mean in case of the ILC5 and HILC5 data, starting with multipoles I > 150. These are most likely caused by 
the extended foregrounds, and point sources remaining in the maps, since we do not apply any sky cuts with these data. We 
over-plot the results for the HILC5 data using the confidence level contours derived from the ILC5 simulations, just to compare 
them with the results obtained for the ILC5 data. 

We also notice that in case of the multipole range- filtered V5 data, (top row in Fig. [2]) some asymmetry is also seen in the 
range £ G [29, 40]. In the same plot it appears that the "northern" hemisphere is anomalous due to the decrement in power as 
compared to the simulations in multipole range £ assy m £ [8, 15]. 

For the case of the multipole, range-filtered results (first and third row in Fig. [2]) we derive the joint "probability of rejecting", 
using a generic multivariate calculus, and for the V5 data we obtain result consistent with the simulations at < 82% confidence 
level for both the maximal hemispherical ratio (^NS,^ max ) and the minimal hemispherical ratio 0~NS,^ max ) regardless of whether 
or not the off-diagonal terms of the covariance matrix are included. This result corresponds to the full range of multipoles up to 
i = 300. 

As a final note to this analysis we point out, that the low significance of the power asymmetry, illustrated in Fig. [2 may result 
from the fact that we analyzed the data in a thin slices through the spherical harmonic space (every two multipoles) up to t = 30. 
Although the full covariance matrix analysis should in principle, be stable to that, it's possible that other binning of the data may 
lead to somewhat different result. In order to check that, we similarly estimated the joint "probability of rejecting" within the 
multipole range £ assy m £ [8,15] for the V5 data, but we found the data to be consistent with simulations at confidence level as 
low as 60%. 

In the following we will constrain the properties of the hemispherical power asymmetry in greater detail. 

IV. MODULATION PARAMETERS ESTIMATION 

In lEriksen et al.l d2007l) an approach for estimating the bipolar modulation parameters was concisely outlined, and was imple- 
mented to obtain the constraints on the modulation parameters. In that work a Gaussian PDF form for all model parameters was 
assumed, except for the modulation orientation axis, for which a flat PDF form was used. As will be shown, the exact shape 
of the likelihood function may have and important effect on the overall significance of the power asymmetry anomaly, so in 
contrast to that work we directly reconstruct the likelihood function using a grid method. 

Since the full parameter space operations are time consuming, we introduce a few assumptions that greatly simplify the 
parameter estimation process. In the next section we discuss them one by one, and provide appropriate justifications. 

A. Bipolar modulation model parametrization 

We generalize the parametrization of the CMB modulation, form the one defined in our previous work dLewll2008h . to account 
for a modulation that is effective only within a requested range of multipoles (f m in, ^max). A CMB observation T(n) of the 
GRF CMB realization (Tcmb(^)) within a bipolar modulation model can be written as: 

T(n) = B(n, fV) * (t C mb(&') (l + M(n')) + F(n')) + N(n)w(h) (2) 
where the modulation field M(n) is defined as: 

M(n) = A^m-n (3) 

where n is a unit vector and M is a bipolar modulation field, oriented along direction rh with amplitude Af™* , which modulates 
the CMB component only within the specified range of multipoles between £ m { n and ^ max . The F(h f ) and N(h f ) denote the 
residual foregrounds and the noise component respectively. The B(n, n') represents the real-space beam convolution kernel of 
the instrument, or any other effective convolution that has been applied to the data. The w(n) is the mask window function 
which can assume either for masked pixels or 1 for unmasked pixels. In case of the ILC5 data w = 1. We will constrain the 
parameters A*™™ and rh in different ranges of multipoles in order to investigate the modulation as a function of scale. 

B. Assumptions 

To facilitate the reconstruction of the multidimensional likelihood function in the following analysis we will rely on three 
assumptions. We assume that (i) the noise in our dataset in the range of multipoles under consideration can be neglected, (ii) that 
the dataset maps are foregrounds free, and (iii) that the residual systematical effects of the modulation-induced change to the 
underlying power spectrum of the CMB does not significantly influences the modulation parameter estimates. In the following 
sections we will discuss each of the assumptions in greater detail. 
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Figure 3: An estimate of the signal-to-noise ratio (SNR) per multipole in the pseudo-power spectrum (C| lgna /Cf olse ) (thick lines) in 
the ILC5 (blue) and V5 (green) dataset and the signal-to-noise ratio of the cumulative variance measurements (thin lines) defined as: 
SNR(^) = E^ 2 ( 2 i + l)C; igna1 / £*, =2 (2?' + l)Q° ise in the ILC5 (blue) and V5 (green) dataset. 



1. Signal to noise ratio 



For the purpose of the analysis we assume that the noise in the CMB dataset (or simulations) has no significant impact on the 
CMB component modulation parameters estimates. The accuracy of this assumption is scale dependent. In Fig. [3] we plot the 
signal to noise ratio (SNR) as deduced from Fig.Q] 

We choose to work with the filtered data at scales below t < 80, where the signal-to-noise ratio per single multipole is 
approximately SNR ~ 10 2 (Fig. [3] thick lines). Note that for the case of a real-space, variance measurements in a map 
composed of a range of multipoles, the SNR is larger due to the cumulative effect (Fig.[3]thin lines). 



2. Foregrounds 

We rely on the foregrounds reduced data of the cleanest WMAP5 channel - V and a conservative sky mask KQ75. As for the 
ILC5 data the residual galactic foregrounds are clearly seen in the map. We drastically reduce their amplitude by a sharp cut 
at the level of ±350/iK (limit deduced from foregrounds free simulations). Of course this doesn't remove the foregrounds but 
somewhat reduces their impact on the regional variance measurements. It should be note d that due to the resid ual foregrounds, 
a caution is recommended by the WMAP team when analysing this map at scales £ > 10 dHinshaw et al .1120071) . We present the 
results of the full sky ILC5 analysis for comparison purposes with the results obtained using an extensively masked V5 data. 
The V5 data should be therefore more reliable in the limit of large multipoles. However it will be interesting to compare the 
results between the two analyses both in the limits of low multipoles, where the ILC map should be reliable, and in the limit of 
large multipoles, where some residual foreground contaminations are present. 



3. Modulation effects to the power spectrum 



The modulation inevitably leads to a change in the underlying power spectrum at all scales, due to the multiplicative dipole 
component. Assuming a modulation orientation along the "z-axis" direction, the modulation field M expressed by a spherical 
harmonic of degree 1 and order is: M(n) = A^ajy^n) = A^^a? cos(0), 2 where a\ = 2 A /|, and the spherical 
harmonic expansion of the modulated map reads: 

T(n) = «r^ m (ft) + 2a/|^ : - Yl a ? cos W y r (*) (4) 

£,m ' £,m 



2 We use the Abramovitz & Stegun (1972) notation convention, and their phase definition of the spherical harmonics throughout the paper. 
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where the first term corresponds to the initial CMB component and the second term corresponds to the modulated component. 
Using the recurrence formula for the associated Legendre polynomials 

{£ - m + VjF^^x) = (2£ + 1)xPP(x) - {£ + m)P^ 1 (x) (5) 

it is straightforward to see that 

cos(0)rr(fi) = ^±l y^ l( n) + i+^y™ l( fi) (6) 

where x = cos(0). 

According to Eq.[6]the modulation leads to redistribution of power of a given multipole £ of the modulation component map 
on to the two neighboring multipoles £ — 1 and t+1 totally removing power from the multipole £. Of course the power in the £i\\ 
multipole is restored by redistribution of power of the £ + 1 and £ — 1 multipoles. In general, the power is redistributed not only 
within the same m mode but also within the mil modes when the modulation orientation is allowed to assume an arbitrary 
orientation. 

It can also be shown that the modulated map has statistically more power than the corresponding non-modulated map. This is 
clearly seen in the modulated map power spectrum (Fig.[4|top panel) as a systematic departure from the initial power spectrum, 
whose magnitude depends on the modulation amplitude A £ £ max . In our analysis we will account for this systematical effect by 
means of a calibration of the map in the real space by the standard deviation (Fig. [4| middle panel). This approach does not 
eliminates the systematical effects, however it reduces them to a sub-percent levels (Fig. |4l bottom panel) for the modulations 
that are cosmologically relevant. 

We study these systematical effects using a sample of 1 000 full sky ILC5 simulations, which we modulate to various extent. 
Then we reconstruct the average modulated power spectra, and compare it to the average non-modulated power spectra (Fig.©. 
Note that e.g. for modulations A < 0.6 the systematical effects of the modulation are smaller than 1% after a proper re-calibration 
of the map. In practice the larger modulations will be much stronger penalized due to the violation of the statistical isotropy, 
than due to the systematical deviations from the initial, fiducial power spectrum. It is also worth noting that the deviation from 
the initial power spectrum in the case of smaller, and much more relevant modulations is negligible, after variance recalibration 
(green curves in the bottom plot of Fig. [4]). 

Since we will be working with selected, filtered ranges of multipoles (in order to investigate the modulation hypothesis as a 
function of scale) the multipoles at the upper and lower ends of that range will have their power significantly suppressed, due 
to the modulation-induced power-leakage outside the chosen range. To improve the effectiveness of the calibration process, 
for a given multipole range of interest [£1^2] we will calibrate the modulation- altered simulations, and the data by standard 
deviations, calculated on the non-modulated, filtered maps, without the outer-most multipoles: i.e. calculated within range 
[£i + 1, £2 — 1]. This improves the statistical consistency by several percent to the case when the calibration is done within the 
full considered multipole range \£\ , £2] . This procedure has been actually used for the results given in Fig ffl 

The bottom panel of the Fig. 0| suggests that the fiducial power spectrum model used by lEriksen et al.l d2007l) : 
^modulated j Qftd _ a (/// Q )& with additional freedom allowed by the tilt parameter in the reference power spectrum could result 
in smaller residuals, however in principle, a possibly large running could also be needed in order to account for the residual, 
systematical discrepancies. 

We also check that the same kind of relation (as depicted in Fig. [4]) with very similar amplitude of the systematical effects is 
obtained for the case of data de-modulation (see. section HVCl) . 

We will test the accuracy of the method to reconstruct the modulation parameters in section [Vj 



C. Method 



Our method is based on measurements of variances in two opposing hemispheres, in a sample of 3 000 GRF V5 and ILC5 
simulations and in the corresponding data. We use 96 different directions in the northern hemisphere, that define the axial- 
symmetry axis of the two hemispherical regions. The directions are defined by the pixel centres of the HEALPIX pixelization 
scheme with the resolution parameter n s = 4. In case of the analysis with V5 data the hemispherical regions are defined outside 
the KQ75 sky mask, and therefore the number of available pixels in these regions (Nk) may vary from one orientation to another. 

In each region k = [1, N r ] (here N r = 2) we measure: 

r k = 4 (7) 
where o\ is the variance of the CMB within /cth region and o 2 is the variance of the whole map (outside the KQ75 sky mask in 
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Figure 4: Top panel: modulation-induced systematical changes to the fiducial power spectrum as inferred from the average power spectrum 
from a sample of 1000 simulations. The overall shift, approximately by a constant factor in all multipoles is clearly seen. Middle panel: average 
power spectra of the modulated maps, calibrated in real space, by the standard deviation. Bottom panel: residual, fractional systematical 
deviation from the original (non-modulated) power spectrum. 



case of the V5 data). For the set of N r regional measurements we define the corresponding x 2 value as: 

N r 

X 2 = E( r * " (rr)) 2 /Var(r s r ) (8) 

k=l 

Note that given that we rely on the local variance measurements it is justified, to the limits to which CMB represents a GRF 
realization, and to the extent to which a cosmic covariance is unimportant, that we neglect the off-diagonal terms in the 
covariance matrix, and treat the regional variance measurements as independent variates. Although we realize that the r/~ 
quantities defined in Eq. [7] should in principle follow a Fischer F-distribution, it is not clear which distribution a sum given in 
Eq. [8] should follow, and as such we reconstruct the likelihood function using a generic y 2 distribution 3 . 

We seek for the best fit between the data and GRF simulations in terms of the hemispherical variance distributions (ri, T2) 
by de-modulating the data under assumptions given in Sect. IIVB1 This approach is therefore a non-standard one, due to the 
fact that generically it's the simulations that are being fit to the data, rather than the data to simulations. However within 
the approximations given in Sect. HVBl it is possible to reverse the process, by demodulating data, while retaining a formal 
correctness and allowing thereby to avoid a time-consuming processing of large number of simulations for each cell of the 



3 We will discuss and take into account possible consequences of this approximation in section PVTC! 
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parameter space. Under no-noise and no-foregrounds assumptions we rewrite the the Eq.[2]as: 

T(n) = B(n,n / )^(T CM B(n / )(l + ^(n / )))^(n) (9) 

It is clear that apart from the beam smoothing effects, the de-modulation of the observed map T(n) is simply a division by the 
factor (l + M(fi / )). In order to account for the beams, using few spherical harmonics transformations (SHT), we pre-process 
the V5 data and simulations as follows: 

1. downgrade the simulations/data to resolution n s = 128 

2. SHT analysis of the full sky V5 maps to obtain a™ V5 coefficients, and deconvolve them with the average V channel beam 
transfer function 

3. SHT synthesis using cl™ V5 nobeam coefficients to obtain a map including the first 128 multipoles (in order to account for 
the power leakage from cut sky in point[4]) 

4. apply KQ75 sky mask to remove the foregrounds in the deconvolved map (there are no foregrounds in the simulations) 

5. SHT analysis of the cut- sky deconvolved map to obtain a™ V5 nobeam cutsky coefficients up to ^ max = 80 

We store the final set of a™v5, nobeam, cutsky coefficients to produce the beam free maps for any requested range of multipoles. 
Note that the two beams of the VI and V2 WMAP channels are practically identical, however we still average between them to 
deconvolve the ILC maps. In the limit of the highest, considered multipoles the operation of deconvolution has an impact of few 
percent, as compared with the convolved power spectrum, as can be inferred from the shape of the beam transfer function. Also 
note, that since we operate in the signal dominated regime there is no danger to artificially blow up the high-£(^ 80) multipoles. 

De-smoothing by the V band beam transfer function, leads to a power increase at t = 80 by about 10%. This can be easily 
estimated from the transfer function itself, since the transmittance for £ = 80 is about 95%, which in the deconvolved power 
spectrum translates onto an increase by a factor of 1/0.95 2 =~ 1.1 

We recall that the noise in the case of V5 is, of course, present in the simulations. The approximation of "no noise" only 
means that we assume that the estimates of the modulation parameter, that modulates the pure CMB component, is not much 
altered by the fact that we are actually deconvolving noisy observations, rather than a pure CMB component, which in general 
does not make sense unless the signal strongly dominates the noise. 

In case of the ILC5 data we create an effective "beam" transfer f unction by dividing the fit to the ILC5 power spectra (see.|Ill 
for details) by the fiducial best fit AC DM model power spectrum dDunkley et al.l d2008l) generated using the mean likelihood 

parameters): i.e. bf = yJcj LCfit /Cf d . We divide the ILC5 af coefficients (t < 80) by this function to match the ILC5 map 
power spectrum to the pure CMB component power spectrum. We will analyze this data with the 3000 GRF signal-only, full 
sky CMB simulations, generated with the best fit fiducial power spectrum Cf d4 . 

As a result, such preprocessed data sets and simulations (apart from the cut sky effects which are identical in the two cases) 
are consistent with the fiducial, best fit, theoretical power spectrum, which we verify experimentally. We can rewrite the Eq. [9] 
as 

^nobeam (n) = T CM B(A)(1 + M(n))w(h) (10) 

On the basis of our preprocessed data set we assured that the inferred modulation parameter will correspond to the modulation 
of the pure CMB component (to the extent where the assumptions given in section HVBl are valid). 

We seek for the best-fit modulation map (1 + M(n)) such that if the observations (T no b e am(n)) are divided by it, the resulting 
map will yield the best consistency with the GRF simulations in terms of the statistics given in Eq.[S] 

D. Parameter space 

As was mentioned in the previous section, we use 96 different directions in the northern sky, that define a set of orientations 
of the N r regions (N r = 2 for hemispherical regions). The set of regions uniformly covers the whole sphere. The 96 directions 
define our search space, and the corresponding search parameter that we will call n s = {1..96}. 



4 Actually for the data preparation process we only generate the GRF aj 1 coefficients and then generate maps for any requested range of multipoles 
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Additionally, we use 192 directions over the full sky, that define the orientation of the de-modulation axis m = {1..192}. The 
directions are defined by the pixel centres of the HEALPIX pixelization scheme of resolution parameter n s = 4. Those that are 
localized in the northern hemisphere overlap with the directions defining the regions orientations. These directions define our 
modulation orientation space. 

In the most general case we probe the likelihood function for the modulation amplitudes in range Af™* G [0.0,0.2] with step 

A = 0.01, and in range A l £°£ G [0.2, 0.3] with step A = 0.02, and in range A l £°* G [0.3, 0.5] with step A = 0.05, and in range 

Aj max G [0.5, 0.7] with step A = 0.1. These values define our modulation amplitude space. 

As will be shown in Sect. [VT] including the large modulations (A > 0.5) mostly explores completely unimportant regions of 
the likelihood function, which is why our grid in this region is much sparser. In general however, the amount of the possible 
hemispherical variance asymmetry in the GRF simulations depends on both: the underlying power spectrum shape, and the 
selected range of multipoles. 

Additionally we perform search in different bins of multipoles (Cin,Cax). The range of the multipoles tested is summarized 
in TableD 



Table I: Summary of the tested multipole bins. Note that as explained in section IIVB 31 in the actual analysis we discard the outermost 
multipoles of the considered multipole ranges. For clarity, in each cell, we explicitly write down the filtered multipole ranges used. The 
numbers in round brackets indicate the percentage of the variance within a considered range out of the total power in the best fit ACDM model, 
calculated as explained in the text. The numbers given in square brackets indicate the percentage of variance within a considered range of 
multipoles out of the total CMB signal variance within the first 80 multipoles of the fiducial power spectrum Cf d . 



2 
6 
19 
29 
39 
59 



3 -6(6.8) [21.2] 3- 
7 



20 

19(13.8) [43.1] 
; 19 (7.0) [21.9] 



3 - 
7- 
20 



30 

29 (16.9) [52.6] 
29(10.1) [31.4] 
- 29 (3.0) [9.5] 



3 - 
7- 
20 
30 



40 

39 (19.4) [60.4] 
39 (12.6) [39.2] 

- 39 (5.6) [17.3] 

- 39 (2.5) [7.8] 



3 - 
7 - 
20 > 
30 
40 



60 

59 (23.8) [74.2] 
59 (17.0) [53.0] 

- 59(10.0) [31.1] 

- 59 (6.9) [21.6] 

- 59 (4.4) [13.8] 



3 - 
7 - 
20 > 
30' 
40 
60 



80 

79 (27.9) [87.2] 
79 (21.1) [66.0] 

- 79(14.1) [44.1] 

- 79(11.1) [34.6] 

- 79 (8.6) [26.8] 

- 79 (4.2) [13.0] 



According to the CMB WMAP5 best fit ACDM model 5 our considered range of the multipoles: i.e. i < 80 make up for 
only about 32% of the total power in this model (of which cumulative variance we calibrate to unity at the maximal computed 
multip ole number of £ = 2000: a(£ min , £ max ) = Yle2T^%( 2e + l ) c f d / E^°2°( 2 ^ + l)Cf d )- However it was shown in lE^wl 
(120081) that the modulation of A « 0.1 extending all the way up to ^ max = 1024 (at which about 96% of the total CMB power 
is used) is excluded at a high confidence level (> 99% CL). 

Throughout the analysis we work on maps of the HEALPIX resolution n s , which depends on the considered range of multi- 
poles so as to yield the condition: n s > £ max /2. 



E. Parameter estimation 



For each direction from our search space (see section HVD1) we reconstruct the likelihood function for each of the the modu- 
lation parameters values = (A, rh) and for each considered multipole range. As a first step we a perform minimization of the 
likelihood over the search parameter n s , in order to select only the measurements that maximize the possible variance distribu- 
tion anomaly. We next derive the corresponding marginal posterior distributions using flat prior probabilities 11(0 1 M) = const. 
at each cell of our parameter space. Therefore, the maximum likelihood inference will lead to the same results as the maximal 
posterior results since according to the Bayes theorem: 

P(0\M,T(&)) oc£(T(&)\M,0)IL(0\M) (11) 

where V(0\M.,T(p)) denotes the posterior distribution, and C(0\Ai,T(h)) denotes the likelihood of the parameters 6 within 
the hypothesized model M defined in Eq.[3l. 

For the estimates on the modulation amplitude, the marginalized, one-dimensional probability distribution is interpolated using 
cubic spline interpolation, before computing the expectancy value, modal value, and confidence ranges. The marginalization over 
the modulation directions is not performed directly on the grid nodes, but rather on an interpolated (on the surface of the sphere, 
for each value of the modulation amplitude independently) posterior. 



5 http://lambda.gsfc.nasa.gov/data/map/dr3/dcp/params/c_l/wmapjcdm_szjens_wmap5_c 
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For the estimates on the modulation orientation the margina lized, two-dimensional probability distribution is interpolated 
using two-dimensional tension splines on sphere dRenkalfl997h . The marginalization over the modulation amplitude is not 
performed directly on the grid nodes, but rather on an interpolated (for each modulation orientation independently) posterior. 
The interpolation is done using a cubic splines and a dense equi-spaced grid. We also tested and compared the "interpolations" 
using spherical harmonics analysis, followed by FWHM = 14° (the approximate size of the search step) Gaussian smoothing 
and synthesis up to ^ max = 30, and found that the results: i.e. the maximum likelihood value orientations and confidence 
contours are reasonably similar, with those where the fitted tension spline surface was used. The usage of cubic interpolations or 
spherical harmonic approximations in principle can lead to oscillations in the PDF that exceed below the zero value especially 
in the tails of the distribution. In case of one dimensional interpolations we circumvent this problem by replacing a cubic 
interpolation with a piece- wise linear interpolation. Although this step will break the continuity of the first derivative of the 
PDF function, we mostly probe the likelihood function dense enough so that these effects are relatively small and unimportant. 
In case of two-dimensional interpolations we find that the potential oscillations, if exist, are small and happen far outside the 
considered confidence levels. In particular such oscillations would lead to artificial generation of multiple isolated contours for a 
given CL which we generally do not observe and consequently do not consider this to be a problem. The most affected artifacts 
of the applied interpolation are observed for cases where the most preferred modulation amplitude is vanishing or is close to 
zero. In such cases, of course, there is no information on the modulation orientation. 

While deriving our results we choose to rely on dense two-dimensional interpolations using tension splines on sphere rather 
than on, somewhat arbitrary, spherical harmonics analysis approach. Depending on the tension parameter the interpolating 
surface approaches the Delaunay triangulation (linear interpolation) solution for large values of the tension parameter, and cubic 
splines solution for zero-tension parameter. The interpolated surface on an equidistant (in galactic latitude and longitude) grid is 
dense enough so that is could easily be projected without holes onto an equal-pixel- area HEALPIX grid to ease the integration 
over the sphere in pixel space. 

V. TESTS OF THE METHOD 

To test the correctness of the code, and the sensitivity of our approach we use a GRF white noise simulations. The choice 
of the white noise helps to generate a GRF simulation in which the cosmic variance effects, leading to accidental, unequal 
hemispherical power distribution, are suppressed, by giving as much power to high-^ modes, as to the low-^ modes. The 
magnitude of the low-£ modes dispersion is oc y^A^/TV^; i.e is larger than that for the high-^ modes, where Ni e uy\ = 21 + 1 
is the number of a™ coefficients at multipoles £, £' respectively, where i > £' . Therefore equalizing power between different 
multipoles helps to better control the experiment: i.e. correctly interpret the results of the tests with the synthetic data given 
some required input modulation parameters. 

In case of realistic GRF simulations based on the CMB power spectrum, the hemispherical power asymmetry allowed within 
the cosmic variance uncertainty, is much larger due to the fact that the lowest multipoles (with the smallest number of modes) 
make up for the main part of the map's total variance, while the higher £ multipoles, even though more numerous, are strongly 
suppressed. 

a. Full sky tests In the following tests there are no effects from any instrumental beams, nor cut skies. We generate a 
white noise realizations in resolution n s = 64 and modulate them with modulation amplitude of A = 0.1 and with modulation 
axis m oriented at (I, b) = (225°, -27°). 

First, we test the correctness of the code by using an analytical proposal for the PDF of the x 2 values. Since we operate on 
white noise, zero mean and unit variance GRF simulations, its statistical properties are well known, and we therefore approximate 
the x 2 value of the Eq.[S]as: 

X 2 = E k= i, 2 (^-(^)) 2 /Var(ai) (12) 

where 

«>=<^U a„c, ^.jjj^j^ (.3, 

are the expectation value of the mean in the sa mple, and the expected variance of the variances in the sample of variates drawn 
from the GRF filed dKenney & Keepin glll95ll) . The a 2 = 1 is the variance of the Gaussian PDF from which the GRF is drawn, 
and Nk is the number of pixels in kth region. 

We find that the statistics correctly reproduced the initial modulation amplitude and orientation well within the "one-sigma" 
confidence level (Fig. \5jp) in all tested cases. Also, the tests show that using the white noise simulations, the method is able 
to reject the hypothesis of statistical isotropy at a very high confidence level (at ^> 4a CL in this case). The precision of the 
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reconstructed, via interpolation, modulation direction orientation is surprisingly good (~ 2°) given a poor resolution of the 
search: - 14° (Fig.0?) 

We also test the reconstruction of the modulation parameters using a set of 1000 GRF simulations from which we derive 
the averages and variances of the regional variance realizations, and proceed according to equations [7] and [8] for N r — 2. We 
notice an increase in the peakedness of the PDF when the Monte-Carlo ("MC") probed estimates are used, as compared to the 
theoretically ("TH") derived estimates. We speculate that the difference might come from the fact that we assumed the field 
variance value a to be unity in all cases. In practice we will always rely on a sample of 3000 MC simulations for estimates of 
local variance distributions. 

Note that the simulation 3 in Fig. [5^, plotted as a peculiarity, found in one of our tests, traces the correct value of the injected 
modulation (A = 0.1) with an accuracy of about 5%. 

b. Cut-sky tests and other subtleties We further test the stability of the method while varying the number of simulations 
used to probe the hemispherical (regional) mean and variance expectations. We check the dependence on the increase of the 
number of simulations from 1000 to 3000. We perform tests with the KQ75 sky cut, and test sensitivity of the method using 
different number of regions: N r = {2, 10}. The regions for N r > 2 are defined as an axial- symmetric patches, equally dividing 
galactic latitude into a symmetrical about the equator (before rotation) regions. In principle the increased number of regions 
could potentially have impact on the accuracy of the method. We use the effective number of degrees of freedom equal to the 
number of unmasked regions, in order to derive the likelihood value. 

The results of the tests are presented in Figs. \5jp-f. We find that the use of the increased number of simulations does not 
significantly influence the estimates of the reconstructed mode values, nor the shape of the marginalized PDF. Also it is clear 
that the increase of the number of regions, used in the statistics, broadens the marginalized PDFs, making thereby the statistic 
less sensitive. 

Furthermore, we see that the accuracy of the method, in case of the cut-sky maps, is generally found at the level of few, up to 
several (in the worst case) percent of the level of the injected modulation (A = 0.1), which is of the same order as the unknown, 
initial (resulting from a random realization) modulation 6 of our white noise maps. 

We note that the selected and presented "simulation 3" is one of the worst cases found in our tests, and as such, we give more 
attention to it in variety of tests summarized in Fig. 0: where the reconstructed distribution exhibits bi-modality. In general 
however, the simulations result in unimodal distributions, like those depicted in Fig. [5^. 

We find that generally, in the presence of the cut sky, the modulation orientation is correctly reconstructed within 50% to 68% 
CL contours (Fig.[5]f), regardless of the number of regions used in the statistics (2 or 10) however in the worst case simulation 
(as in case of the "simulation 3" with the two-region statistics) it is found as far as within the 95% CL contour. 

We also checked the difference between different statistical approaches: i.e. between maximization 7 and marginalization over 
the modulation orientation. We find that both - the full sky and cut sky tests yield similar, or almost identical results. We will 
show that the situation will not be the same in case of the real CMB data or CMB simulations due to the effects we mentioned 
in the beginning of this section. 

Finally we tested the statistical biases of the method under the cut sky conditions (Fig. \5]f (left) red lines) and found that, within 
the obtained accuracy, no significant statistical bias is noticed. 

c. Summary We find that with our chosen search resolution, the method traces the correct solution to within a few percent 
accuracy for the full sky measurements with respect to the injected modulation amplitude value, and from few up to several 
percent accuracy (< 18%) for the cut sky case, with about 68% of the estimates yielding an accuracy better than ~ 7% (Fig. [2). 
In terms of the absolute errors of the reconstructed value of the modulation amplitude parameter A, for the injected amplitude of 
A = 0.1 the errors are roughly an order of magnitude smaller: ~ 0.005 and ~ 0.007 for the full and cut sky cases respectively. 
In case of the larger number of regions the sensitivity of the method is worsened (eg. for the case of N r = 10) and therefore in 
the following analysis we will only rely on the two-region statistics. 

As for the reconstructed modulation direction, we find that mostly the correct direction is reconstructed within ~ 50% CL 
limits for the two regions statistics, and well within 50% CL limit for 10 regions case for the cut sky and full sky cases. 

It is important to note that even with the white-noise simulations the initial, unknown modulation, resulting from random, and 
unequal distribution of power in the sky is at level of A < 0.005, which is of the same order of magnitude as the accuracy which 
we obtain in the full sky tests. 

As for the reconstructed modulation orientation, mostly the correct direction is found to be within the 50% to 68% CL limits 
in case of the cut sky reconstructions, while the typical angular size of the 50% and 95% CL contours are ~ 20° and ~ 35° 



6 We refer to the initial unknown, accidentally unequal power distribution in a GRF white noise simulation as a "modulation" since it's the modulation amplitude 
that we measure, but of course there's no reason to believe that any modulation effect, as defined in this paper, exists in the GRF simulations. 

7 In case of maximization over the modulation direction orientation, we have used the modal values, found in the fitted, two-dimensional maps of the likelihood 
function, for each modulation amplitude. We found this method to improve the smoothness of the resulting PDF, since our parameter search space is very 
sparse - only 192 directions over the entire sky. 
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Full sky white-noise maps tests with initial modulation: A = 0.1 (Z, b) — (225° , —27°) 



a) marginalized modulation amplitude PDF: 



sim 1: TH, mode={0.1, 0.104} 
sim 2: TH, mode={0.1, 0.0993} 
sim 2: MC, mode={0.1, 0.0999} 
sim 3: MC, mode={0.095, 0.0954} 



0.06 



0.08 



0.10 0.12 0.14 0.16 
modulation amDlitude (A) 



0.18 



0.20 




b) marginalized 50%, 68% and 95% CL modualtion orientation limits 

mode: abl«=^24 r 16] 



0.22 



Cut sky (KQ75) white-noise maps tests with initial modulation: A = 0.1 (/, b) = (225°, -27° 



c) marginalized modulation amplitude PDF: 



sim 3: 2 regs., TH, full sky, mode={0.095, 0.0972} 
sim 3: 2 reg., MC (1000 sim.)., full sky, mode={0.095, 0.0954} 
sim 3: 2 reg., MC (1000 sim.), cut sky, mode={0.12, 0.1175} 
sim 3: 10 reg., MC (1000 sim.), cut sky, mode={0.11, 0.1094} 
sim 3: 10 reg., MC (3000 sim.), cut sky, mode={0.11, 0.1102} 



0.10 0.11 
modulation amplitude (A) 




d) marginalized 50%, 68% and 95% CL modualtion orientation limits 

mode: ab?»(2l&.-20] 




Two-region statistic tests with cut sky (KQ75) white-noise maps and with initial modulation: A = 0.1 (/, b) = (225°, —27°) 
e) marginalized modulation amplitude PDF: f) statistics of the reconstructed ML modulation amplitude values 



sim 4: 2 regs., MC(1000 sim.), cut sky, mode={0.09, 0.0902} 
sim 5: 2 reg., MC (1000 sim.), cut sky, mode={0.105, 0.1054} 
sim 4: 10 reg., MC (1000 sim.), cut sky, mode={0.085, 0.086} 
sim 5: 10 reg., MC (1000 sim.), cut sky, mode={0.095, 0.0932} 
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Figure 5: White noise maps tests of the accuracy of the method to reconstruct the injected modulation amplitude and its orientation, (a) 
Constraints on the modulation amplitude, from three, chosen, full sky, white-noise, modulated maps. The three colors code three selected 
simulations. The dashed lines for each simulation connect points at which the likelihood function was calculated, and the solid lines represent 
a cubic polynomial (spline) interpolation, combined with a linear interpolation in the tails of the distribution so as to avoid the oscillations into 
negative values. The abbreviations "TH" (for the black and the blue-spiky curves) and "MC" (for the green and the broad-blue lines) in the 
legend indicate respectively: the theoretically derived, and Monte-Carlo probed values of the expected means and variances of the regional 
(hemispherical) variances (Eq. [T3j. In the "MC" case a sample of 1 000 simulations is used. The mode values of the probed and interpolated 
PDFs are also given, (b) Limits of the modulation orientation, from one of the modulated, white noise, full sky, GRF realization. The 50% 
(dark blue) , 68% (light blue) and 95% (red) confidence level ranges are plotted. The small red dot indicates the location of the m direction, 
and in the top of the plot the ML modulation direction is given. We obtain similar results in all tested cases, (c) As in panel (a) but for the 
KQ75 cut-sky, modulated, white-noise maps tests. We focus here only at the peculiar, "worst-case" - the "simulation 3", subject to different 
statistical approaches as indicated in the legend. For comparison, the "TH" PDF, and the PDF for the full sky case for this simulation, are 
also replotted. (d) Limits of the modulation orientation for the sky cut "simulation 3", using N r = 10 region statistic. The details are as 
in panel (b). (e) Constraints on the modulation amplitude, from two another (less-extreme than "simulation 3") simulations, derived using 
the two-region statistics (black curves), and the corresponding PDFs for 10 regions statistics (red curves), (f) Left: reconstructed modulation 
PDF mode values versus expectation values from a 100 cut-sky simulations, with 68% CL error bars derived individually from each PDF 
independently, plotted with the corresponding histograms (right). The CL ranges were integrated from the mode value. The red lines indicate 
the mean value of the scatter for each direction, and the shaded area encompasses one standard deviation of the distribution. 
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respectively (Fig. [5]). The full- sky reconstructions CL contours are slightly smaller. Note that the modulation direction is 
reconstructed via interpolation to within a few degrees accuracy for the full sky case with the search resolution of about 14°, 
which is surprisingly good. For the cut sky case the accuracy is approximately at the level of several degrees. 



VI. RESULTS 



A. Modulation amplitude 



The marginalized over the modulation direction, probability distributions (posteriors) of the modulation amplitude parameter 
as derived from the WMAP V5 and ILC5 data are plotted in figures @and[7]respectively for different ranges of filtered multipoles 
(see table J]). We remind the reader that only the analysis involving the V5 data was performed using the KQ75 sky mask. 

TablelUsummarizes the modal values of the distributions as parameter estimates and the 68% and 95% confidence level limits. 



Table II: Constraints on the modulation amplitude parameter for the V5 and ILC5 data. The table contain the modal parameter values and the 
corresponding 68% and 95% (in brackets) confidence level limits. 



2 



V5 data 

7 20 30 40 60 80 

(0.00) 0.00 < 0.11 < 0.29 (0.55) (0.01) 0.12 < 0.23 < 0.29 (0.34) (0.00) 0.07 < 0.19 < 0.24 (0.28) (0.00) 0.05 < 0.16 < 0.21 (0.25) (0.00) 0.00 < 0.00 < 0.15 (0.20) (0.00) 0.00 < 0.12 < 0.15 (0.18) 

(0.07) 0.14 < 0.21 < 0.26 (0.31) (0.00) 0.00 < 0.14 < 0.18 (0.20) (0.00) 0.06 < 0.13 < 0.16 (0.18) (0.01) 0.05 < 0.10 < 0.12 (0.14) (0.02) 0.05 < 0.08 < 0.11 (0.13) 

(0.00) 0.00 < 0.01 < 0.06 (0.12) (0.01) 0.05 < 0.09 < 0.12 (0.15) (0.01) 0.04 < 0.07 < 0.09 (0.11) (0.02) 0.04 < 0.05 < 0.07 (0.09) 

(0.01) 0.06 < 0.10 < 0.14 (0.19) (0.00) 0.00 < 0.05 < 0.07 (0.09) (0.00) 0.00 < 0.00 < 0.04 (0.06) 

(0.00) 0.00 < 0.02 < 0.03 (0.07) (0.00) 0.00 < 0.00 < 0.03 (0.06) 
(0.00) 0.01 < 0.03 < 0.05 (0.07) 

ILC5 data 



(0.00) 0.00 < 0.00 < 0.26 (0.52) (0.07) 0.16 < 0.28 < 0.33 (0.41) (0.06) 0.14 < 0.23 < 0.28 (0.33) (0.06) 0.13 < 0.22 < 0.26 (0.30) (0.03) 0.09 < 0.17 < 0.20 (0.24) (0.03) 0.08 < 0.15 < 0.18 (0.20) 
(0.01) 0.12 < 0.22 < 0.28 (0.31) (0.04) 0.09 < 0.15 < 0.20 (0.24) (0.06) 0.10 < 0.13 < 0.17 (0.20) (0.04) 0.07 < 0.09 < 0.12 (0.15) (0.03) 0.05 < 0.08 < 0.10 (0.12) 
(0.02) 0.06 < 0.10 < 0.13 (0.17) (0.02) 0.05 < 0.07 < 0.09 (0.12) (0.00) 0.02 < 0.04 < 0.06 (0.07) (0.01) 0.02 < 0.03 < 0.04 (0.06) 
(0.00) 0.04 < 0.07 < 0.10 (0.12) (0.00) 0.00 < 0.00 < 0.03 (0.05) (0.00) 0.00 < 0.00 < 0.02 (0.04) 
(0.00) 0.01 < 0.02 < 0.03 (0.05) (0.00) 0.00 < 0.00 < 0.02 (0.03) 
(0.00) 0.00 < 0.00 < 0.03 (0.06) 



In Table HITl we concisely summarize the results of the modulation significance analysis: ie. the analysis in which we derive the 
minimal confidence levels, at which the modulation value of A = cannot be excluded. We specify the expectancy values, mode 
values of the distributions, and the corresponding significance. We choose to calculate the confidence intervals - or rather, since 
we're working on posterior probability distributions, in the nomenclature of the Bayesian language, the credibility intervals, by 
integrating from the modal value, rather than from the expectancy value. 



Table III: Results of the modulation amplitude parameter estimation for the V5 and ILC5 dataset. The table contain the minimal confidence 
levels (as percents) at which the parameter value of A = cannot be excluded (bold face numbers) and the expectancy (in round brackets) and 
the modal (in square brackets) values of the corresponding distributions. See also Fig. [8] 



2 

6 
19 
29 
39 
59 

tiin\-^n 

2 

6 
19 
29 
39 
59 



7 20 
57.5 (0.23) [0.11] 96.1 (0.20) [0.23] 
99.5 (0.19) [0.21] 



V5 data 
30 



40 



90.9 (0.15) [0.19] 88.7 (0.13) [0.16] 
90.2 (0.11) [0.14] 94.5 (0.11) [0.13] 
25.8 (0.05) [0.01] 96.8 (0.08) [0.09] 
97.1 (0.10) [0.10] 



60 

4.6 (0.10) [0.00] 
97.0 (0.08) [0.10] 
97.2 (0.06) [0.07] 

86.0 (0.05) [0.05] 

61.1 (0.03) [0.02] 



ILC5 data 

7 20 30 40 

0.4 (0.20) [0.00] 99.1 (0.24) [0.28] 99.1 (0.20) [0.23] 99.3 (0.19) [0.22] 
96.6 (0.19) [0.22] 99.2 (0.14) [0.15] 99.9 (0.13) [0.13] 
97.9 (0.09) [0.10] 99.2 (0.07) [0.07] 
95.7 (0.07) [0.07] 



60 

98.7 (0.14) [0.17] 
99.7 (0.09) [0.09] 
95.6 (0.04) [0.04] 
0.4 (0.02) [0.00] 
78.1 (0.02) [0.02] 



80 

86.4 (0.09) [0.12] 
99.4 (0.08) [0.08] 

99.7 (0.05) [0.05] 
0.2 (0.03) [0.00] 
0.3 (0.02) [0.00] 

88.8 (0.04) [0.03] 

80 

98.8 (0.12) [0.15] 

99.6 (0.08) [0.08] 

98.7 (0.03) [0.03] 
0.4 (0.02) [0.00] 
0.5 (0.01) [0.00] 
0.2 (0.03) [0.00] 



To visualize these results, we plot the estimated modal values of the posterior distributions as a function of considered £ m i n 
and ^ max values (Fig. [8]). For each multipole bin, in Fig. [8] we also indicate the minimal confidence level (see Table (TTTJ) it takes 
to keep the modulation A = i.e. the non-modulated, isotropic model, as an viable option. 
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Figure 6: Constraints on modulation- amplitude parameters as a function of multipole bin considered, from the V5 data. From the bottom to the 
top we plot the data reducing filtering of the large scale multipoles. For each PDF the corresponding expectancy value is marked by vertical 
dashed lines for bins £ e [3, 6], £ G [3, 19], £ £ [3, 29], and dash-dotted line for bins £ e [3, 39], £ e [3, 59] and £ e [3, 79]. Within each group 
the increasing line width corresponds to increasing value of £ ma ^- Only every 100th point of the interpolated, marginalized PDF was plotted, 
so the data points do not correspond to the actual grid nodes. We have truncated the plot at A=0.4 to maximally expose the most interesting 
regions, while keeping the same scale throughout all panels. 
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A 

Figure 7: Constraints on modulation- amplitude parameter as a function of filtered multipole bin from the ILC5 data. Details as in Fig. [6] 
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Figure 8: Results of the modulation amplitude parameter estimation from the V5 data (top panel) and ILC5 data (bottom panel). Modal values 
of the posterior distributions are plotted for different bins of multipoles filtered out from the dataset. The error bars represent the confidence 
limits that include the A = value. The corresponding confidence levels (given in percents) are given as annotations by the data points. 



It is clear that the modulation amplitude depends on the considered multipole range and hence on scale. Is is generally 
seen that large values of the best-fit modulations mostly come from the large scales, while for high multipole bins, the best-fit 
modulations are much smaller. As explained in section [VI within GRF realizations this is somewhat expected due to the nature 
of cosmic variance effects. 

Looking at distributions in Figs. [6] and [3 it is apparent that the modulation for i G [7, 19] (and also t G [3, 19]), in the V5 
data, (blue curves (crosses) in Fig. [6]) is strongest and most significant, as it takes the confidence level as high as ~ 99.5% to 
include the A = value (see.Table [Till and Fig. [8]). The appearance of some asymmetry in this range also seems consistent with 
the results presented in Section [III1 The range of multipoles i G [20, 29] of the V5 data does not seem to prefer any modulation 
value as its modal value is almost zero, and consequently while increasing values of 4nax, for £ m i n = 7, the overall significance 
falls, as this multipoles bin is included, but then systematically increases as larger multipoles are added, which is consistent with 
the shape of the PDF functions preferring some non-zero modulation for higher multipole bins like £ G [30, 39]. 

The best-fit modulations for i G [20, 39] and t G [30, 39] range from A = 0.07 to A = 0.10 and exclude the isotropic model, 
depending on the data, at confidence level of about 96% to 99%, while the best-fit modulation for t G [7, 39] in the ILC5 data, 
with the modal amplitude of A = 0.13, exclude the isotropic model, at confidence level as high as 99.9%. We will further test 
the significance of these results in section IvFc! 



18 



Note also that, some of the marginalized PDFs exhibit bi-modality (eg. the PDF corresponding to the multipole range 
£ G [3, 59] - magenta-line (V) in Fig. [6]). This bi-modality results from the marginalization itself, and is not observed in the 
full non-marginalized distribution. Since the likelihood function does not depend on the orientation of the modulation axis for 
the modulation amplitude A = 0, while it does depend on the modulation orientation very strongly for modulations A ^> 0, 
the likelihood surface shall tend to be peaky for large values of A, and flat for A « 0. Hence, depending on how strong the 
preference of some direction happens to be, it is possible to accumulate in the marginalization process a second peak (the second 
mode) out of somewhat less-preferred, but constant at certain level, likelihood values along A « direction. As a result the 
aforementioned range £ G [3, 59] of V5 data yields a small significance in Table [TTT] (in terms of rejecting an isotropic model). 
We have also processed the results by using maximization over the modulation orientations instead of marginalizations, and 
as expected, the maximized PDFs are unimodal and more strongly exclude the non-modulated, isotropic models, however we 
choose the more conservative, and more correct method of marginalizing over the non plotted dimensions. 

From Fig. [8] it is easy to see that the modulation amplitude estimates are mostly similar between the two datasets, and that 
generally the amplitude of the modulation decreases with increasing multipole number £ in the two datasets. A common feature 
between the two datasets is that for high multipole bins £ G [40, 59], £ G [40, 79] and £ G [60, 79] the best-fit modulation 
amplitudes are small or zero, despite that the amount of variance carried by bin eg. £ G [40, 79] is as high as ~ 27% of the 
total variance carried by the full range of considered multipoles (i.e. from £ = 3 to £ = 79, see. table J]). The multipoles range 
£ G [29, 40] instead participate to the total variance only by ~ 8%, and consequently, the the best-fit modulation of this range, 
estimated to be A = 0.10 (A = 0.07) for the V5 (ILC5) data, is effectively destroyed, as the higher multipole bins are included, 
most likely due to simply dominating power in the added multipole bins. 

There are few significant differences between the datasets as well. Firstly, we notice, that the ILC5 estimates are generally 
slightly, but systematically larger from the V5 estimates. Also, in particular, the full sky ILC5 posteriors, for multipole bins 
£ G [7, 29] and £ G [20, 29] strongly prefer some non- vanishing modulation amplitude, in contrast to the V5 data (compare 
cyan-diamonds in Fig. [6] and This results in an almost constant significance of excluding A = as higher multipole bins are 
being included (see first three rows of the Table [TTT] in section for ILC5 data). 

It is interesting to note a small difference in range £ G [3, 6] (see table [ill]) in which the ILC5 slightly prefer a vanishing 
best-fit modulation. Note that among our considered multipole bins, the ILC5 should be reliable basically only within this 
lowest range. In contrast, some non-zero preferred modulation is obtained with the V5 data; however the value is still largely 
consistent with the vanishing modulation at confidence level as low as ~ 58%. 

Some differences between the datasets are of course expected due to the cut-sky effects, which preclude filtering of exactly 
the same range of multipoles due to the power leakage effects in case of V5 data. Also caution is needed in the interpretation 
of the ILC5 data for higher multipole bins, as residual foregrounds in the regions around the Galactic center, may have some 
impact on the results. In particular, these residual foregrounds might be responsible for the significant alteration of the shape of 
the PDF function in the multipole bin £ G [20, 29] (cyan-diamonds in figures[6]and U} towards an increased significance in favor 
of non-isotropic models. 

B. Modulation orientation 

We now focus on the modulation orientation as a function of our chosen multipole bins, as specified in table J] The maximum 
likelihood modulation orientations are summarized in table HVl 

The best-fit orientations found in the analysis depend in general on the considered range of the multipoles. In particular, we 
see that the hemispherical power asymmetry, as measured here by the modulation orientations, generally tend to shift from larger 
galactic latitudes to the smaller galactic latitudes, as dat a of higher multipole bins are processed or cumulatively added up as 
tax value increases. This was previously also noticed in Eewl(l2008h but using a different method. This effect seems to be seen 
in the first and second rows (as moving from the left to the right) for each of the data in table [IVJ or by comparing vertically in 
columns the first three rows as the low multipoles are removed from the analysis. This is also seen directly in the distributions 
plotted in Fig. [9l where it can be easily deduced in which bins the effects of small or vanishing modulations can be ignored. 
However we note that the effect is not present in every bin, and at most of order of few tens of degrees, and therefore we refrain 
from making any far-reaching speculations based on it. 

While analyzing the distributions in Fig.|9]note that, those of them that correspond to very small or vanishing modal values of 
modulation amplitudes, have a very extended confidence level contours that cover large fraction of the sphere. These naturally 
result from a very flat likelihood function, and therefore any inference based on these cases is largely speculative and irrelevant. 
This is however expected, because for a vanishing modulation amplitude, there is no information on its orientation either. 

We note that it is possible there is some degree of correlation between the plots for a given row, resulting from an cumulative 
effect of adding higher multipole bins. The possible changes to the resulting distribution will jointly depend on the modulation 
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Figure 9: Constraints on the modulation orientation from the V5 and ILC5 data. For each considered multipole bin we plot the confidence 
regions corresponding to the 50% (green), 68% (yellow), and 95% (red) confidence levels, based on the interpolated maps of the posterior 
distributions. In each map the maximum likelihood orientation and the corresponding dipole plane are indicated using a light blue dot and line 
respectively. Additionally, for comparison, the ecliptic south pole and ecliptic plane are plotted in dark blue. In the top of each panel we give 
the galactic coordinates of the maximum likelihood solution. The arrangement of the panels is consistent with the cells in tables[TJ|llJ andHTIl 
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amplitude and orientation in the added bin, but also on the amount of variance carried by the that bin as specified in table U 

In table|lV]we plot the directions, in galactic coordinates, of the maximum posterior values found in the modulation orientation 
analysis. Out of curiosity we also provide the angular separation of these directions from the ecliptic south pole to check for any 
possible extra alignments. 



Table IV: Results of the modulation orientation parameter estimation for the V5 and ILC5 data. The table contains the galactic coordinates 
of the maximum posterior modulation orientation and (in square brackets) the relative angular distance to the south ecliptic pole at (/, b) = 

(276.4°, -29.8°). 



i-mm \<-max 

2 



V5 data 

7 20 30 40 60 

(281°, -20°), [11°] (233°, -54°), [39°] (234°, -46°), [37°] (225°, -47°), [43°] (220°, -42°), [46°] 
(278°, -68°), [38°] (242°, -55°), [35°] (213°, -42°), [52°] (225°, -35°), [43°] 
(187°, 2°), [90°] (205°, -8°), [70°] (213°, 0°), [67°] 
(168°, -19°), [96°] (193°, 1°), [85°] 
(202°, 8°), [81°] 

ILC5 data 

7 20 30 40 60 

(205°, -67°), [56°] (226°, -53°), [43°] (224°, -56°), [45°] (225°, -55°), [44°] (226°, -56°), [44°] 
(281°, -57°), [27°] (264°, -48°), [21°] (256°, -43°), [21°] (257°, -42°), [20°] 
(236°, 19°), [62°] (231°, -6°), [49°] (241°, 2°), [46°] 
(191°, -18°), [77°] (225°, 1°), [58°] 
(311°, -34°), [30°] 
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-34°), [45°] 
-31°), [45°] 
,1°),[64°] 
-2°), [67°] 
, 1°), [49°] 
12°), [43°] 

80 

-53°), [43°] 
-42°), [19°] 
-5°), [42°] 
-15°), [49°] 
-18°), [17°] 
-47°), [46°] 



C. Modulation significance 

In the previous section we have shown that, for some multipole ranges, the reconstructed, marginalized probability distribution 
function of the modulation amplitude, excludes the vanishing modulation value (A = 0) at a very high confidence level. It is 
important to ask whether this result is really robust, and whether or not we should reject the standard isotropic model of the 
Universe, at least, at some of the scales: i.e. those corresponding to the distributions with the strongest modulation detections, 
and the highest non-zero modulation significances. In particular, at least three ranges are of most concern: £ £ [7, 19], and 
£ e [7, 79] for which the A = can be excluded at 99.5% and 99.4% CL respectively, using the KQ75 sky-cut V5 data, and 
where the modulation parameters are constrained to be within ranges (0.07)0.14 < 0.21 < 0.26(0.31) at 68% (95%) CL and 
(0.02)0.05 < 0.08 < 0.11(0.13) at 68% (95%) CL respectively. Also the aforementioned range £ e [7, 39], for which the A = 
can be excluded at 99.9% CL using the full- sky ILC5 data, and where the modulation parameter is constrained to be within 
range (0.06)0.10 < 0.13 < 0.17(0.20) at 68% (95%) CL. 

What we have done in the previous sections, is that we have estimated the best-fit modulation parameters (amplitude and 
orientation) with respect to the average from large amount of GRF simulations. Using average from large number of simulations 
ensures that we compare the data to really isotropic distributions, as any deviations from the statistical isotropy, even those 
resulting from the cosmic variance, will be averaged out. Although the measurements also quantified the allowed magnitude 
of deviation from the ideal isotropy, allowed within the cosmic variance, via the standard deviation in the \ 2 tests, there were 
number of explicit, or implicit assumptions or simplifications made on the way, like for example the usage of the y 2 distributions, 
or ne glecting the cosmic covariance effects, or residual foregrounds to name few. 

In iLewl d2008|) we have performed the full covariance matrix analysis in two hemispherical regions with the same search 
parameter space as detailed in Section HVD1 We used 1000 GRF simulations (500 for covariance matrix estimation and 500 for 
probing the PDF of the underlying y 2 distribution) and another 1000 simulations, modulated with amplitude of A = 0.114. The 
simulations were filtered up to £ max = 40. We found that on average about 8% of the GRF simulations exhibited a more unusual 
power distributions, than those found in the modulated simulations. In this paper, although we use the same density in the search 
space, we have improved somewhat the method by using smooth interpolations. 

In order to further test the robustness, and the significance of the power asymmetry anomalies, and circumvent all possible 
imperfections of the method, in the following we will pursue a similar test. We process 100 GRF simulations of the V5 data 
through our parameter estimation pipeline, and compare the results with the real data. Such approach should always be an 
ultimate test of the robustness, as it must give correct results independently from the assumptions taken in the method. 

We will focus on the aforementioned multipole bins: £ e [7, 19] and £ G [7, 79] of the V5 data, and £ e [7, 39] of the ILC5 
data. Within these ranges the power asymmetry seems to be very strong and very significant (see. table[III]and figures and [7]). 
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Figure 10: In the left-hand side panels we plot a comparison of the reconstructed, marginalized distributions of the modulation amplitudes 
from 100 simulations of the WMAP V5 data (top and middle panels) and WMAP ILC5 data (bottom panel). The WMAP data are plotted 
using red, thick lines (triangles). Only every 100th point of the interpolated, marginalized PDF is plotted. In the right-hand side panels we plot 
the corresponding histograms of the expectancy values and modal values derived from these distributions. The WMAP data values are marked 
with vertical lines. The plotted ranges yield the strongest, and most significant, plausible hemispherical anomalies in the data as inferred from 
the analysis in section IvTa! 
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In Fig.[K)|we plot the results of the modulation parameter estimation for all tested simulations along with the WMAP data. 
While it is clear that most of the simulations do not prefer any significant (if any at all) modulation amplitude values, at least 
few simulations, in our sample, yield modulations that are stronger in the considered range of multipoles, than those found in 
the data. Also, from the shape of the PDF it is easy to infer that the significance of rejecting A = in these few cases will be 
even larger than in the case of the selected, most anomalous results from the WMAP data. 

We find that 7 out of 100 V5 simulations yield stronger best-fit (modal) modulation amplitudes, and 6 of them also yield a 
more significant rejection of the A = parameter value, than the V5 data in the range £ G [7, 19]. Similarly, for the range 
£ G [7, 79] 5 simulations yield stronger and more significant best-fit modulation values. 

Consequently, we conclude, that the significance, as inferred simply from integrating the PDF (as given in tableHTTJ) is not quite 
robust. In light of these results we estimate the overall significance of possible modulation signals in the analyzed WMAP CMB 
maps at the level of about ~ 94% t o ~ 95% de pending on the particular range of multipoles. This remains greatly consistent 
with our previous results reported in lLewl(l2008l) forthe same data. 

As for the ILC5 data we find that three simulations out of 100 yield a more significant rejection of A = hypothesis, and 
curiously, only one in those three also yields a stronger modulation within the multipole range £ G [7, 39] . Therefore, the 
corresponding overall significance of the power asymmetry, in this particular multipole range, is still as high as about ~ 99%. 
We note however, that since we did not use any sky masks in this case, this result can probably be safely considered as somewhat 
overestimated, as any residual galactic foregrounds are likely only to increase the level of the hemispherical power asymmetry, 
rather decrease it. 



VII. DISCUSSION 

The results outlined in sections [VI] indicate that different multipole ranges yield a different best fit modulation value, and that 
the modulation orientation also slightly varies from one multipole range to another. 

In particular, the best-fit modulation orientation dependence when higher multipole bins are included cumulatively is not as 
strong, as when the added multipole bins are considered individually. Generally, we notice that within the best-fit orientations, 
that also yield a large modulation values (A > 0.1) the high-£ multipole bins prefer a close galactic plane orientation, while the 
low-£ multipole bins rather prefer orientations with larger galactic latitudes. 

The analysis of the modulation amplitude within few multipole bins yielded a large, best-fit modulation amplitudes, that seem 
to significantly reject the isotropic Universe model (with A = 0). However as much as few in one hundred GRF simulations, 
processed as data, also yielded a similar or larger modulation values, and also excluded the A = hypothesis at yet even higher 
confidence levels, than in the case of the V5 data. This effectively reduces the overall significance at which the isotropic model 
of the Universe can be rejected, down to only about 94% or ~ 95% using the V5 data in the range £ G [7, 19], and £ G [7, 79] 
respectively. 

We therefore pursued the analysis of the modulation signals in a two partially complementary ways. While the first approach 
addresses the question of "how large and how significant is the best-fit modulation of the data?", the second approach quantifies 
"how consistent is the best fit-modulation as compared with the GRF simulation expectations?". The second approach should be 
more robust since it is free of any, possibly inaccurate, assumptions that could result in underestimation of the size of the errors 
in the statistic, and in the result lead to spurious detections. These problems are effectively eliminated in a direct comparison 
with the GRF simulations. 

Curiously the ILC 5 data in the multipo le range £ G [7, 39] still seem to be anomalous at a high CL of about 99%; level almost 
as high as quoted in dEriksen et alJl2007h . However contrary to that work, we have not applied any sky masks to this data, and 
therefore these results, given here only for comparison purposes with the V5 data, should still be treated with caution. 

It would be interesting to perform similar analysis using the ILC5 data but with included sky cut, and to check the dependence 
of the analysis while varying the sky cuts from less to more aggressive. Also, it could be interesting to check the robustness 
and the significance in other multipole ranges than those two, tested in section PvTc! In principle, it would also be interesting 
to include other available renditions of the ILC maps, to see the stability of the modulation to different foregrounds cleaning 
pipelines. We defer these issues for possible future work. 

VIII. CONCLUSIONS 

We performed tests of the hemispherical power asymmetry found in the CMB WMAP data for different bins of multipoles in 
two ways. 

First, we introduced a statistic that searches for the orientation of two opposing, hemispherical regions that maximize or 
minimize the hemispherical variance ratio, and compared these with the expectations from the GRF simulations. We found that 
the maximal asymmetry revealed this way is found within a multipole range £ G [8, 15], with the southern hemisphere having 
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larger variance than the northern hemisphere. When these results are compared to the GRF simulations, the northern hemisphere 
appears to be suppressed below the average expectation. 

Secondly, we have introduced and tested a new metho d for measuring the power asymmetry in the CMB data, as quantified 
within a bipolar modulation model (IGordon et al ] l2005b . For the first time we constrained the modulation parameters as a 
function of various multipole bins. For each multipole range, we obtained the constraints on the the modulation amplitude and 
orientation. On the basis of the data sets, analyzed up to the maximal multipole £ max = 80, i.e. the WMAP five-year inverse 
noise co-added, KQ75 sky cut map from the V channel (V5), and the five-year, full-sky, foregrounds cut ILC map (ILC5) we 
found that: 

(i) generally the modulation amplitude decays as higher multipole bins are cumulatively added or independently analyzed, 

(ii) the best fit modulation amplitude is small A < 0.03 and insignificant for multipoles beyond £ w 40 

(iii) the most anomalous signals in terms of the modulation amplitude and its significance come from multipole range t G [7,19], 
and i G [7, 39] in the V5 and ILC5 data respectively. For these ranges the significances of rejecting the isotropic cosmological 
model are 99.5% and 99.9% respectively and the constraints on the best fit, (PDF modal) modulation values are: (0.07)0.14 < 
0.21 < 0.26(0.031) and (0.06)0.10 < 0.13 < 0.17(0.20) at 68% (95%) CL respectively. 

Focusing on the two selected multipole ranges we performed an additional tests of the significance using GRF simulations 
processed as data, and found that similar or stronger and more significant (in terms of rejecting the isotropic model) modulation 
values are obtained in 6 (1) cases in 100 simulations, which decreases the overall significance of the power asymmetry in the 
CMB down to 94% (99%) in V5 (ILC5) data respectively. To complement the results in the limit of high multipoles as well, we 
additionally tested the range t G [7, 79] of the V5 data that also yields a strong and significant (99.4%) best-fit modulation value 
- (0.02)0.05 < 0.08 < 0.11(0.13) at 68% (95%) CL - but when this result was compared with the GRF simulations the effective 
significance is again decreased down to about 95%. 

Although the significance in case of the ILC5 data is still rather high, we warn that the results in this case were obtained 
without any sky cut, and therefore the asymmetry significance can be overestimated due to residual foregrounds. 

Finally we note that a further analysis of the significance in terms of comparison with GRF simulations of other multipole 
ranges would be interesting, as well as analysis of the power asymmetry in the ILC data as a function of different sky cuts. 
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